/*
Copyright (c) 2009 David Bucciarelli (davibu@interfree.it)

Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:

The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/

#ifndef _GEOMFUNC_H
#define	_GEOMFUNC_H

#include "geom.h"
#include "simplernd.h"

#ifndef SMALLPT_GPU
float path_medio = 0.0f;
int l_path_medio = 1;
int currentSamplePath = 0;


/*La funzione torna 0 se la sfera i-esima non viene colpita, altrimenti t*/
static float SphereIntersect(
#ifdef GPU_KERNEL
OCL_CONSTANT_BUFFER
#endif
	const Sphere *s,
	const Ray *r) { /* returns distance, 0 if nohit */
	Vec op; /* Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0   OVVERO L'EQUAZIONE DELLA SFERA, scritta in forma vettoriale*/
	vsub(op, s->p, r->o); //s->p è il centro della sfera

	float b = vdot(op, r->d);
	float det = b * b - vdot(op, op) + s->rad * s->rad;
	if (det < 0.f) //Viene risolta come se fosse una equazione di secondo grado, quindi se il determinante è <0 non abbiamo soluzione,
		           //in senso ottico il raggio non colpisce la sfera
		return 0.f;
	else
		det = sqrt(det);

	float t = b - det;
	//Se tutti e due gli t (+ o -) sono negativi vuol dire che la sfera è dietro al raggio, quindi non si vede
	if (t >  EPSILON)
		return t;
	else {
		t = b + det;

		if (t >  EPSILON)
			return t;
		else
			return 0.f;
	}
}
//In questo modo si crea un sampling dell'emisfero (?)
static void UniformSampleSphere(const float u1, const float u2, Vec *v) {
	const float zz = 1.f - 2.f * u1;
	const float r = sqrt(max(0.f, 1.f - zz * zz));
	const float phi = 2.f * FLOAT_PI * u2;
	const float xx = r * cos(phi);
	const float yy = r * sin(phi);

	vinit(*v, xx, yy, zz);
}

/*Funzione che verifica se il raggio r intercetta una sfera della scena */
static int Intersect(
#ifdef GPU_KERNEL
OCL_CONSTANT_BUFFER
#endif
	const Sphere *spheres,
	const unsigned int sphereCount,
	const Ray *r,
	float *t,
	unsigned int *id) {
	float inf = (*t) = 1e20f; //Punto lontano, intersezione a distanza t

	unsigned int i = sphereCount;
	for (; i--;) { //Controlla ogni sfera della scena, una per volta, salvando l'id dell'intersezione più vicina
		const float d = SphereIntersect(&spheres[i], r);
		if ((d != 0.f) && (d < *t)) {
			*t = d;
			*id = i;
		}
	}

	return (*t < inf);
}

//Funzione che controlla se il punto colpito di una sfera non colpisca altre sfere
static int IntersectP(
#ifdef GPU_KERNEL
OCL_CONSTANT_BUFFER
#endif
	const Sphere *spheres,
	const unsigned int sphereCount,
	const Ray *r,
	const float maxt) {
	unsigned int i = sphereCount;
	for (; i--;) {
		const float d = SphereIntersect(&spheres[i], r);
		if ((d != 0.f) && (d < maxt))
			return 1;
	}

	return 0;
}

/* La funzione di occupa di vedere per ogni fonte luminosa
 * cerca lo shadow ray che colpisce tale luce dal punto colpito
 * e non interseca nessuna altra sfera lungo il percorso
 */

static void SampleLights(
#ifdef GPU_KERNEL
OCL_CONSTANT_BUFFER
#endif
	const Sphere *spheres,
	const unsigned int sphereCount,
	unsigned int *seed0, unsigned int *seed1,
	const Vec *hitPoint,
	const Vec *normal,
	Vec *result) {
	vclr(*result);
	float wo=0;
	int visible=0;

	/* For each light */
	unsigned int i;
	for (i = 0; i < sphereCount; i++) {
#ifdef GPU_KERNEL
OCL_CONSTANT_BUFFER
#endif
		const Sphere *light = &spheres[i];
		if (!viszero(light->e)) {
			/* It is a light source */
			Ray shadowRay;
			shadowRay.o = *hitPoint;

		   do{
			    wo = 0; visible=0;
			    /* Choose a point over the light source */
			    Vec unitSpherePoint;
			    UniformSampleSphere(GetRandom(seed0, seed1), GetRandom(seed0, seed1), &unitSpherePoint);
			    Vec spherePoint;
			    vsmul(spherePoint, light->rad, unitSpherePoint);
			    vadd(spherePoint, spherePoint, light->p);

			    /* Build the shadow ray direction */
			    vsub(shadowRay.d, spherePoint, *hitPoint);
			    const float len = sqrt(vdot(shadowRay.d, shadowRay.d));
			    vsmul(shadowRay.d, 1.f / len, shadowRay.d);

			    wo = vdot(shadowRay.d, unitSpherePoint);
			    if (wo < 0.f) {
			    	visible = 1;
			    	wo = -wo;
				    // continue;
			    } //else /* It is on the other half of the sphere */


			    /* Check if the light is visible */
			    if(visible){
			        const float wi = vdot(shadowRay.d, *normal);
			        if ((wi > 0.f) && (!IntersectP(spheres, sphereCount, &shadowRay, len - EPSILON))) { //Quindi non interseca nessuna altra sfera lungo il percorso
				        Vec c; vassign(c, light->e);
				        const float s = (4.f * FLOAT_PI * light->rad * light->rad) * wi * wo / (len *len);
				        vsmul(c, s, c);
				        vadd(*result, *result, c);
			        }
			    }
		   //printf("wo=%f\n",wo);

		   }while(!visible);
		}
	}
}


static void SampleLightsOld(
#ifdef GPU_KERNEL
OCL_CONSTANT_BUFFER
#endif
	const Sphere *spheres,
	const unsigned int sphereCount,
	unsigned int *seed0, unsigned int *seed1,
	const Vec *hitPoint,
	const Vec *normal,
	Vec *result) {
	vclr(*result);

	/* For each light */
	unsigned int i;
	for (i = 0; i < sphereCount; i++) {
#ifdef GPU_KERNEL
OCL_CONSTANT_BUFFER
#endif
		const Sphere *light = &spheres[i];
		if (!viszero(light->e)) {
			/* It is a light source */
			Ray shadowRay;
			shadowRay.o = *hitPoint;

			/* Choose a point over the light source */
			Vec unitSpherePoint;
			UniformSampleSphere(GetRandom(seed0, seed1), GetRandom(seed0, seed1), &unitSpherePoint);
			Vec spherePoint;
			vsmul(spherePoint, light->rad, unitSpherePoint);
			vadd(spherePoint, spherePoint, light->p);

			/* Build the shadow ray direction */
			vsub(shadowRay.d, spherePoint, *hitPoint);
			const float len = sqrt(vdot(shadowRay.d, shadowRay.d));
			vsmul(shadowRay.d, 1.f / len, shadowRay.d);

			float wo = vdot(shadowRay.d, unitSpherePoint);
			if (wo > 0.f) {
				/* It is on the other half of the sphere */
				continue;
			} else
				wo = -wo;

			/* Check if the light is visible */
			const float wi = vdot(shadowRay.d, *normal);
			if ((wi > 0.f) && (!IntersectP(spheres, sphereCount, &shadowRay, len - EPSILON))) {
				Vec c; vassign(c, light->e);
				const float s = (4.f * FLOAT_PI * light->rad * light->rad) * wi * wo / (len *len);
				vsmul(c, s, c);
				vadd(*result, *result, c);
			}
		}
	}
}



/* --RadiancePathTracing--
 * E' la funzione che si occupa di trovare le intersezioni tra il raggio in input e i vari oggetti della scena(sfere)
 * Una volta trovato un punto colpito dal raggio si deve capire il tipo di materiale dell'oggetto colpito,
 * in modo da calcolare i raggi di riflessione e rifrazione.
 * */
static void RadiancePathTracing(
#ifdef GPU_KERNEL
OCL_CONSTANT_BUFFER
#endif
	const Sphere *spheres,
	const unsigned int sphereCount,
	const Ray *startRay,
	unsigned int *seed0, unsigned int *seed1,
	Vec *result) {

	Ray currentRay; rassign(currentRay, *startRay);
	Vec rad; vinit(rad, 0.f, 0.f, 0.f);
	Vec throughput; vinit(throughput, 1.f, 1.f, 1.f);
    double p=0;
    float rand=0;
    float path_tmp =0;

	unsigned int depth = 0;
	int specularBounce = 1;

	//New_media=((old_media*old_contatore)+depth_attuale)/new_contatore), dove new_contatore=old_contatore+1).

	//Indica il livello di profondità del path tracing
	for (;; ++depth) {
		// Removed Russian Roulette in order to improve execution on SIMT
//		if (depth > 6) { //Se l'algoritmo non ha raggiunto almeno 6 passi di profondità non viene fermato
//			*result = rad;
//			return;
//		}

		float t; /* distance to intersection */
		unsigned int id = 0; /* id of intersected object */
		if (!Intersect(spheres, sphereCount, &currentRay, &t, &id)) {
			*result = rad; /* nel caso in cui non incontri nessun oggetto lungo il percorso, ritorna */
			path_tmp = path_medio;
			if(currentSamplePath == 0)
			{
			    path_medio = depth; currentSamplePath =1;
			}
			else
			{
	        	path_tmp = path_medio;
			    path_medio = ((path_medio*l_path_medio)+depth)/(l_path_medio+1);
			   // printf("Sto eseguendo (%.3f * %d)+%d/%d e ottengo %.3f\n",path_tmp,l_path_medio,depth,l_path_medio+1,path_medio);sleep(1);
			    l_path_medio++;
			   // printf("La lunghezza del path medio: %d\n",l_path_medio);sleep(1);
			 }
			return;
		}

#ifdef GPU_KERNEL
OCL_CONSTANT_BUFFER
#endif
		const Sphere *obj = &spheres[id]; /* the hit object, trovato con Intersect */

       Vec f=obj->c;
       if (f.x > f.y && f.x > f.z){
    	   p = f.x;
       }
       else if(f.y > f.z){
    	   p = f.y;
       }
       else p = f.z;
       // p = (f.x + f.y + f.z)/3;

       //SPARA!!!! RR
       rand = GetRandom(seed0,seed1);
       if(rand > p && depth>6){
    	   *result = rad;
    	   //printf("Sono uscito con %f>%f mentre la profondita era:%d\n",rand,p,depth);
    	   *result = rad; /* nel caso in cui non incontri nessun oggetto lungo il percorso, ritorna */
           if(currentSamplePath == 0) {
        	   path_medio = depth; currentSamplePath =1;
           }
           else
           {
        	   path_tmp = path_medio;
        	   path_medio = ((path_medio*l_path_medio)+depth)/(l_path_medio+1);
        	  // printf("Sto eseguendo (%.3f * %d)+%d/%d e ottengo %.3f\n",path_tmp,l_path_medio,depth,l_path_medio+1,path_medio);sleep(1);
        	   l_path_medio++;
        	  // printf("La lunghezza del path medio: %d\n",l_path_medio);sleep(1);
           }

    	   return;

       }
//       else
//       {
//    	   f.x=f.x*(1/p);
//    	   f.y=f.y*(1/p);
//    	   f.z=f.z*(1/p);
//
//       }


        //Calcola punto colpito(intersecato) e normale
		Vec hitPoint;
		vsmul(hitPoint, t, currentRay.d);
		vadd(hitPoint, currentRay.o, hitPoint);

		Vec normal;
		vsub(normal, hitPoint, obj->p);
		vnorm(normal);

		//Quanto un punto viene colpito dal raggio, dobbiamo capire il tipo di materiale per calcolare il raggio di rifrazione (DOT)
		const float dp = vdot(normal, currentRay.d);

		Vec nl;
		// SIMT optimization
		const float invSignDP = -1.f * sign(dp);
		vsmul(nl, invSignDP, normal);

		/* Add emitted light */
		Vec eCol; vassign(eCol, obj->e);
		if (!viszero(eCol)) {
			if (specularBounce) {
				vsmul(eCol, fabs(dp), eCol);
				vmul(eCol, throughput, eCol);
				vadd(rad, rad, eCol);
			}

			*result = rad;
			path_tmp = path_medio;
			if(currentSamplePath == 0)
			{
				path_medio = depth; currentSamplePath =1;
			}
			else
			{
				path_tmp = path_medio;
				path_medio = ((path_medio*l_path_medio)+depth)/(l_path_medio+1);
				//printf("Sto eseguendo (%.3f * %d)+%d/%d e ottengo %.3f\n",path_tmp,l_path_medio,depth,l_path_medio+1,path_medio);sleep(1);
				l_path_medio++;
				//printf("La lunghezza del path medio: %d\n",l_path_medio);sleep(1);
			 }
			return;
		}



		if (obj->refl == DIFF) { /* Ideal DIFFUSE reflection */
			specularBounce = 0;
			vmul(throughput, throughput, obj->c);

			/* Direct lighting component */

			Vec Ld;
			//Per un tipo di riflessione diffusa si deve calcolare la luce in tutte le direzioni
			SampleLightsOld(spheres, sphereCount, seed0, seed1, &hitPoint, &nl, &Ld);
			vmul(Ld, throughput, Ld);
			vadd(rad, rad, Ld);

			/* Diffuse component */

			float r1 = 2.f * FLOAT_PI * GetRandom(seed0, seed1);
			float r2 = GetRandom(seed0, seed1);
			float r2s = sqrt(r2);

			Vec w; vassign(w, nl);

			Vec u, a;
			if (fabs(w.x) > .1f) {
				vinit(a, 0.f, 1.f, 0.f);
			} else {
				vinit(a, 1.f, 0.f, 0.f);
			}
			vxcross(u, a, w);
			vnorm(u);

			Vec v;
			vxcross(v, w, u);

			Vec newDir;
			vsmul(u, cos(r1) * r2s, u);
			vsmul(v, sin(r1) * r2s, v);
			vadd(newDir, u, v);
			vsmul(w, sqrt(1 - r2), w);
			vadd(newDir, newDir, w);

			currentRay.o = hitPoint;
			currentRay.d = newDir;
			continue;
		} else if (obj->refl == SPEC) { /* Ideal SPECULAR reflection */
			specularBounce = 1;

			//In questo caso tramite il dot product otteniamo il raggio speculare
			Vec newDir;
			vsmul(newDir,  2.f * vdot(normal, currentRay.d), normal);
			vsub(newDir, currentRay.d, newDir);

			vmul(throughput, throughput, obj->c);

			rinit(currentRay, hitPoint, newDir);
			continue;
		} else {
			specularBounce = 1;

			Vec newDir;
			vsmul(newDir,  2.f * vdot(normal, currentRay.d), normal);
			vsub(newDir, currentRay.d, newDir);

			Ray reflRay; rinit(reflRay, hitPoint, newDir); /* Ideal dielectric REFRACTION */
			int into = (vdot(normal, nl) > 0); /* Ray from outside going in? */

			float nc = 1.f;
			float nt = 1.5f;
			float nnt = into ? nc / nt : nt / nc;
			float ddn = vdot(currentRay.d, nl);
			float cos2t = 1.f - nnt * nnt * (1.f - ddn * ddn);

			if (cos2t < 0.f)  { /* Total internal reflection */
				vmul(throughput, throughput, obj->c);

				rassign(currentRay, reflRay);
				continue;
			}

			//Se l'oggetto non ha riflesso completamente il raggio al suo interno
			//si calcola riflessione o rifrazione

			//Si calcola il raggio rifratto
			float kk = (into ? 1 : -1) * (ddn * nnt + sqrt(cos2t));
			Vec nkk;
			vsmul(nkk, kk, normal);

			Vec transDir;
			vsmul(transDir, nnt, currentRay.d);
			vsub(transDir, transDir, nkk);
			vnorm(transDir);

			//Riflessione o rifrazione usanto il termine di Fresnel
			float a = nt - nc;
			float b = nt + nc;
			float R0 = a * a / (b * b);
			float c = 1 - (into ? -ddn : vdot(transDir, normal));

			float Re = R0 + (1 - R0) * c * c * c * c*c;
			float Tr = 1.f - Re;
			float P = .25f + .5f * Re;
			float RP = Re / P;
			float TP = Tr / (1.f - P);

			if (GetRandom(seed0, seed1) < P) { /* R.R. */
				vsmul(throughput, RP, throughput);
				vmul(throughput, throughput, obj->c);

				rassign(currentRay, reflRay);
				continue;
			} else {
				vsmul(throughput, TP, throughput);
				vmul(throughput, throughput, obj->c);

				rinit(currentRay, hitPoint, transDir);
				continue;
			}
		}
  }

}

static void RadianceDirectLighting(
#ifdef GPU_KERNEL
OCL_CONSTANT_BUFFER
#endif
	const Sphere *spheres,
	const unsigned int sphereCount,
	const Ray *startRay,
	unsigned int *seed0, unsigned int *seed1,
	Vec *result) {
	Ray currentRay; rassign(currentRay, *startRay);
	Vec rad; vinit(rad, 0.f, 0.f, 0.f);
	Vec throughput; vinit(throughput, 1.f, 1.f, 1.f);

	unsigned int depth = 0;
	int specularBounce = 1;
	for (;; ++depth) {
		// Removed Russian Roulette in order to improve execution on SIMT
		if (depth > 6) {
			*result = rad;
			return;
		}

		float t; /* distance to intersection */
		unsigned int id = 0; /* id of intersected object */
		if (!Intersect(spheres, sphereCount, &currentRay, &t, &id)) {
			*result = rad; /* if miss, return */
			return;
		}

#ifdef GPU_KERNEL
OCL_CONSTANT_BUFFER
#endif
		const Sphere *obj = &spheres[id]; /* the hit object */

		Vec hitPoint;
		vsmul(hitPoint, t, currentRay.d);
		vadd(hitPoint, currentRay.o, hitPoint);

		Vec normal;
		vsub(normal, hitPoint, obj->p);
		vnorm(normal);

		const float dp = vdot(normal, currentRay.d);

		Vec nl;
		// SIMT optimization
		const float invSignDP = -1.f * sign(dp);
		vsmul(nl, invSignDP, normal);

		/* Add emitted light */
		Vec eCol; vassign(eCol, obj->e);
		if (!viszero(eCol)) {
			if (specularBounce) {
				vsmul(eCol, fabs(dp), eCol);
				vmul(eCol, throughput, eCol);
				vadd(rad, rad, eCol);
			}

			*result = rad;
			return;
		}

		if (obj->refl == DIFF) { /* Ideal DIFFUSE reflection */
			specularBounce = 0;
			vmul(throughput, throughput, obj->c);

			/* Direct lighting component */

			Vec Ld;
			SampleLights(spheres, sphereCount, seed0, seed1, &hitPoint, &nl, &Ld);
			vmul(Ld, throughput, Ld);
			vadd(rad, rad, Ld);

			*result = rad;
			return;
		} else if (obj->refl == SPEC) { /* Ideal SPECULAR reflection */
			specularBounce = 1;

			Vec newDir;
			vsmul(newDir,  2.f * vdot(normal, currentRay.d), normal);
			vsub(newDir, currentRay.d, newDir);

			vmul(throughput, throughput, obj->c);

			rinit(currentRay, hitPoint, newDir);
			continue;
		} else {
			specularBounce = 1;

			Vec newDir;
			vsmul(newDir,  2.f * vdot(normal, currentRay.d), normal);
			vsub(newDir, currentRay.d, newDir);

			Ray reflRay; rinit(reflRay, hitPoint, newDir); /* Ideal dielectric REFRACTION */
			int into = (vdot(normal, nl) > 0); /* Ray from outside going in? */

			float nc = 1.f;
			float nt = 1.5f;
			float nnt = into ? nc / nt : nt / nc;
			float ddn = vdot(currentRay.d, nl);
			float cos2t = 1.f - nnt * nnt * (1.f - ddn * ddn);

			if (cos2t < 0.f)  { /* Total internal reflection */
				vmul(throughput, throughput, obj->c);

				rassign(currentRay, reflRay);
				continue;
			}

			float kk = (into ? 1 : -1) * (ddn * nnt + sqrt(cos2t));
			Vec nkk;
			vsmul(nkk, kk, normal);
			Vec transDir;
			vsmul(transDir, nnt, currentRay.d);
			vsub(transDir, transDir, nkk);
			vnorm(transDir);

			float a = nt - nc;
			float b = nt + nc;
			float R0 = a * a / (b * b);
			float c = 1 - (into ? -ddn : vdot(transDir, normal));

			float Re = R0 + (1 - R0) * c * c * c * c*c;
			float Tr = 1.f - Re;
			float P = .25f + .5f * Re;
			float RP = Re / P;
			float TP = Tr / (1.f - P);

			if (GetRandom(seed0, seed1) < P) { /* R.R. */
				vsmul(throughput, RP, throughput);
				vmul(throughput, throughput, obj->c);

				rassign(currentRay, reflRay);
				continue;
			} else {
				vsmul(throughput, TP, throughput);
				vmul(throughput, throughput, obj->c);

				rinit(currentRay, hitPoint, transDir);
				continue;
			}
		}
	}
}

#endif

#endif	/* _GEOMFUNC_H */

